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ABSTRACT 


The Segmented Mirror Telescope (SMT) housed at the Naval Postgraduate School is a unique, 
state-of-the-art optical instrument built to explore new technologies needed for future space- 
based telescopes. A discrete Fourier transform wavefront reconstruction technique developed 
at Lawrence Livermore National Laboratory is discussed in this thesis as applied to a hexagonal 
aperture. A Fourier domain implementation of a spatial-frequency modal controller for a simple 
spring-mass model of a deformable mirror surface is provided by this thesis. This technique 
avoids more difficult time-domain solutions, is computationally efficient and scalable to much 


larger multi-input, multi-output systems than the SMT. 
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Executive Summary 





The Segmented Mirror Telescope (SMT) was built to develop new techniques for rapid con- 
struction of optical space-based telescopes. This system has been transferred to the Naval 
Postgraduate School (NPS) and is now hosted by the Spacecraft Research and Design Cen- 
ter (SRDC). As a part of this initiative, the SRDC now is the Adaptive Optics Center of Ex- 
cellence for National Security and is actively pursuing research focused on the SMT and other 
optics systems. 


The most complex problem of optical systems is constructing large primary mirrors of sufficient 
optical quality. Traditional monolithic mirrors can take up to four years to construct and test. If 
the mirror is segmented and split into several smaller mirrors, the process becomes much more 
rapid. The mirror segments for the SMT were constructed in approximately 18 months. These 
smaller mirrors can also be combined to form much larger mirrors than could have been built 


by the traditional monolithic processes. 


The trade-off in this new technique is an increased complexity of the system. The monolithic 
mirror is built to act as a coherent reflector across the entire mirror surface. The larger size 
also means that it is not as susceptible to damage during the rocket launch into orbit. For the 
segmented system, the mirrors must be aligned correctly. If this process is not done properly and 
the mirror height between segments varies too much, the resulting imaging is worse. Phasing 
of the telescope, the highly sensitive alignment process, must be done to ensure this does not 


occur. 


In addition to the mirror surface being segmented, the rear surface also has surface-parallel 
piezoelectric actuators. These actuators can correct for mirror imperfections. This arrangement 
is referred as an active-optics deformable mirror system. This allows for increased variance 
tolerances in the mirror manufacturing process since the mirror can be corrected when the ac- 
tuators are energized. This technique lowers costs associated with construction and testing but 


adds complexity to the control system. 


In contrast to the active-optics system, SMT is an adaptive-optics system. The telescope in- 
cludes a Shack-Hartmann Wavefront Sensor (WES), which provides information about the qual- 
ity of the light reaching the camera. The light may lose phase coherency and become less planar. 
This occurs from either atmospheric disturbances, structural jitter and/or mirror imperfections. 


Each source of error is summed to provide the total error. The WES provides the feedback 


Xlll 


mechanism for the mirror actuators. A closed control loop would be suited for this situation 
to minimize wavefront error by correcting the light to have a planar wavefront at the camera 


detector. 


In order to correct the light, an estimate of the wavefront must be formed. An implementation 
of the Poyneer algorithm is developed in MATLAB to handle other aperture shapes, such as 
the hexagons used in the SMT. The sampling geometry for the Shack-Hartmann is applied 
to the algorithm. The algorithm provides the information necessary to be used in a real-time 
implementation of the wavefront reconstruction such that it could be used in a feedback control 


loop. 


A novel controller design that uses the two-dimensional fast Fourier transform (FFT2) is pro- 
posed in this thesis to solve this problem. First, a simple model is developed where a deformable 
mirror surface is treated as a grid of masses and springs. The actuators have surface normal 
forces applied to individual masses. The equation to describe this arrangement is a triple convo- 
lution in both the space and time domains. Because of the complexity of the coupled equations, 
the FFT2 is employed. In the frequency domain, the problem transforms into a much simpler 


arrangement of many uncoupled, independent equations. 


After developing the mirror model, we designed the controller model. Instead of the traditional 
design of one controller for each actuator, one controller for each spatial mode of actuators 
is used. Each controller then provides actuator commands for every actuator. The inverse 


transform is performed to recover the specific commands for each actuator. 


The mirror model and the controller are implemented in this thesis using MATLAB Simulink. 
The simulation takes advantage of the powerful discrete filter blocksets to efficiently process 
the transfer functions in parallel. The mirror aperture size is set to be 256x256 to show the 
algorithm is scalable to much larger apertures than are currently used by state-of-the-art optical 


systems. 


A closed-loop feedback control simulation is performed with time-varying phase wavefront 
information. The simulation is run to verify the design of the controller and the wavefront error 
decrease as the actuators affect the mirror surface. The results are presented in the thesis. This 
form of a controller can be used on the actual SMT after determining the actual coefficients of 
the mechanical structure of the system. 


The implementation of a discrete Fourier transform (DFT) wavefront reconstruction and de- 


XIV 


formable mirror model and controller is developed in this thesis. These tools will be useful 
in the future research of the SMT and can be applied to an actual system for experimental re- 
sults. Any techniques that improve the imaging quality of the SMT can be applied to other 
telescope systems to ensure this technology continues to progress with the goal of providing 


higher quality image products for less cost. 
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CHAPTER 1: 


Introduction 





Traditional space-based optical telescopes employ a single, large primary mirror structure. This 
mirror has to be painstakingly manufactured to minimize deformation errors across the surface. 
The testing may take months to years, depending on the program requirements. These large 
mirrors are costly to manufacture, test and schedule for spacecraft. The astronomy community 
has developed a number of technologies that allow for segmented mirrors to replace single, large 
mirrors. This technology is now being developed by the space-based telescope community for 


the next generation of telescopes [1]. 


The Naval Postgraduate School (NPS) hosts a Segmented Mirror Telescope (SMT), originally 
produced by the National Reconnaissance Office (NRO) as a proof-of-concept experimental 
spacecraft telescope [2]. It incorporates many aggressive concepts into a single working unit to 


test and learn about the manufacturing of such an instrument. 


The primary mirror is segmented into six equal hexagons arranged in the same manner as a 
traditional mirror. The segments are paired with each other and have two sets of hinges. In this 
arrangement, the segmented mirror can be folded up for launch and mechanically unfold and 
effloresce the segments. This allows for a smaller launch configuration, which increases options 


for launch vehicles. 


In addition to the primary mirror being segmented, it also has actuators across the rear of the 


Figure 1.1: SMT segmented mirror orientation. 


surface. There are three types of actuators: coarse, fine and face sheet. These actuators can 
modify the mirror surface in such a way as to correct for abnormalities in the approaching light. 
The coarse and fine actuators are used as a calibration to give the primary mirror a high-quality 
optical surface. These actuators make the system an active-optics system. With a wavefront 
sensor these distortions are measured, and a control system adjusts the face-sheet actuators to 
compensate. In this manner, the telescope can also be considered an adaptive-optics system, 
which is an improvement over active optics. What makes this unique is that the primary mirror 
has actuators, whereas most adaptive-optics systems make use of a smaller deformable mirror 


such as a fold mirror that is after the primary in the optical path [3]. 


The National Security Space Strategy [4] outlines the importance of systems like the SMT. 
Similar systems offer “unprecedented advantages in national decision-making, military oper- 
ations and homeland security.’ The United States must maintain the benefits afforded by our 
exquisite systems in the evolving strategic environment. Developing technology that improves 
the design of the SMT directly supports these goals. 


1.1 Purpose 

Some of the background fundamental concepts required for a space-based optical telescope that 
functions as an adaptive-optics system are developed in this thesis. The SMT has a unique 
geometry that requires finesse in application of these concepts. The concepts are introduced 
for a simple square mirror and then applied to the SMT and modeled. Simulation results are 
presented. The improvements made to the SMT model and controller will result in more accu- 
rate analysis in future space telescope research. This work benefits government agencies and 


companies interested in space telescopes. 


The purpose of this research is to improve the quality of the imaging results of the telescope. The 
ideal wavefront is planar at the camera sensor; although, the wavefront can be distorted from 
the ideal geometry for a number of reasons. Atmospheric mixing and motion can cause non- 
uniform delays in the wavefront, a common problem for ground-based astronomical telescopes. 
For spacecraft telescopes, the structural jitter from the motion of changing the satellite pointing 
direction can deform the optics. For Earth-observation satellites, this is a great concern because 
the pointing direction is changed often during operation. In this research, a sample wavefront 
is used in the discussion. The three-dimensional image of the wavefront is shown in Figure 1.2. 


The red areas indicate the leading wavefront, and the blue areas indicate the lagging wavefront. 





Figure 1.2: Projection of the sample wavefront. 


1.2 Overview 

The background theory of wavefront reconstruction, which is a process employed to correct the 
optical system, is provided in Chapter 2. This process is applied to a square and circular aper- 
tures, which was developed in previous work [5,6]. Wavefront reconstruction is then discussed 


for the hexagonal aperture and the unique characteristics of the SMT. 


The two-dimensional fast Fourier transform (FFT2) is used to model the SMT in Chapter 3. 
Use of the FFT2 leads to a computationally efficient algorithm that yields good results. This 
technique yields decoupled dynamics and a more efficient control design. This improves the 
scalability of the processing such that larger apertures can be computed without large increases 


in computational time. 


The implementation and simulation work done for this thesis is covered in Chapter 4. Wavefront 
reconstruction from example generated data is presented. The technique can also be applied to 


noisy slope data from actual systems. The Deformable Mirror (DM) controller using the discrete 


Fourier transform (DFT) is also described, and the technique of building the actual model and 
its simulation are also presented. The transfer function equivalents are considered to reduce the 


order of the system. Sample phase data is used to illustrate the performance of the system. 


The summary of the work and conclusions with recommended follow on research for the SMT 


are discussed in Chapter 5. 


1.3. Research Contributions 

Square and circular apertures are discussed in existing literature and analyzed in great detail 
[5-8]. Although upcoming telescopes use hexagonal apertures, they are not widely treated 
in the literature. The required modifications for the wavefront reconstruction algorithm of a 


circular aperture to correctly model a hexagonal aperture are provided. 


Zernike polynomials are used commonly in literature and practice to characterize systems and 
wavefronts [9]. Although the polynomials are particularly suitable to model optical wavefronts, 
they are computationally inefficient. An alternate technique that is computationally efficient 
is desired for realtime, closed-loop operations of a DM. A spatial frequency wavefront recon- 


struction technique is developed as a viable alternate to Zernike polynomials. 


A further contribution is also made in the control system design. In a traditional control design, 
commands are generated for each actuator individually in the system. This does not scale 
efficiently as the number of actuators increases. A modal controller is developed to correct for 


wavefront distortions in the spatial frequency domain. 





CHAPTER 2: 


Wavefront Reconstruction 





2.1 Theory 


Wavefront reconstruction is the mathematical process of using the sensed wavefront gradient 
information and determining an estimated equivalent wavefront. An approaching wavefront has 
very fine perturbations in the phase which can be observed in the resulting images as distortions. 


The differences in phase can be specified by using the optical path length [10] 


Sp ise (2.1) 
ny 

where 6¢@ is the difference in phase, A is the wavelength of light, n is the index of refraction 
of the media and éd is the minute difference in the physical path length, expressed as a lump 
sum of all of the distortions. The quantity nd is the optical path length. In the laboratory 
environment, the medium is air. The optical path length is usually considered for the center 
wavelength of interest, in this case within the visible light band. Thus, the correction cannot be 
applied uniformly for all wavelengths, as each has a different optical path length. 


These differences reduce the resolution and clarity on the resulting image. To correct for the 
phase differences, the mirror is properly reformed by a set of actuators. This can be seen in 
Figure 1.2, where the mirror needs to provide a negative phase adjustment for the red areas by 
pushing the mirror back and a positive phase adjustment for the blue areas by pushing the mirror 
forward, so that when the light arrives at the mirror surface, the wavefront reflection is planar. 


The end result is improved imaging. 


In the actual implementation of Equation (2.1), the phase @ is an estimate that cannot fully 
eliminate all errors. The phase estimates @ (x,y) are computed across a two-dimensional grid of 
values on the basis of observed local gradients. The gradients are sampled at a low spatial fre- 
quency across the aperture and, as a consequence of the sampling, can only reconstruct specific 
low-order modes of the wavefront. These gradients are then processed to determine a wave- 
front that has local gradients of the same form as the original wave. The device that samples 
the wavefront is called the wavefront sensor [11]. There are several different wavefront sensors; 


the one that is discussed in this thesis is the Shack-Hartmann sensor [12], which is a grid of 





Figure 2.1: Square grid of phase points that contribute to gradients in (a) Hudgin geometry 
and (b) Fried geometry. Shack-Hartmann sensors are shown in red dashed lines with hexagon 
lenslets. 


small lenses that are side-by-side. Each lenslet focuses the received light to a point on the focal 
plane of a secondary camera. If the point does not appear at the center of the lenslet optical axis 
on the camera focal plane, then the light has wavefront variation. The position of these points 


relative to their expected centered position yields the information on the phase gradient. 


There are two common sensor geometries developed by Hudgin [7] and Fried [8] shown in 
Figure 2.1. These geometries vary by the location of the Wavefront Sensor (WFS) relative 
to the phase points, which are abstract locations in the optical plane that have the equivalent 
gradients between adjacent phase points as registered by the sensors. The Hudgin geometry 
is easier to work with initially in simulation but has the drawback that the sensor’s lenslets 
have overlap, which is not commonly done in actual implementation. Fried geometry is more 
commonly used to model the actual Shack-Hartmann sensors. The SMT uses hexagon-shaped 
lenslets for its sensors. These pack together densely without overlap and focus all available 
light on the sensor detection focal plane. 


Wavefronts can be reconstructed using a variety of basis functions. Some reconstructions use 
Zernike polynomials, which are particularly suitable to describe the common optical character- 
izations of astigmatism, coma, defocus and others [9]. Although a full expansion is computa- 


tionally inefficient, only a few terms are needed to characterize a wavefront. 


An alternate decomposition by the DFT yields the desirable property of being computationally 
efficient and still only providing a few dominant coefficients [6]. The work in this thesis uses 


this technique. 


2.2 Square Aperture 


The square aperture is a primitive geometry that allows for a simplified reconstruction algo- 
rithm. The design assumes that the actuators and phase points are co-located. Thus, the actu- 
ators directly influence the phase points (and, thereby, the gradients). For a NxN grid of phase 
points [m,n], Hudgin’s geometry model defines the gradients as the difference of the neighbor 


phase points: 


Sx[m,n] = o|m+ 1,n] — o[m,n] (2.2a) 
sy[m,n] = o[m,n+ 1] — o[m,n]. (2.2b) 


Alternatively, Fried’s model defines the gradients as the average of the neighbor phase points: 


s,{m,n] = $(olm-+ 1,n] — o[m,n| + 6[m+ 1,n+4+ 1] — o[m,n+ 1}) (2.3a) 


Sy|m,n] = $(Olmn-+1] — [m,n] + o[m+ 1,n+ 1]—o[m+1,n)). (2.3b) 


In either case, the result is two sets of slope data that are an (N — 1)x(N — 1) grid of values, 
rather than NxN because the slopes can only be computed between data points. Slopes are 
defined as the local gradients separated into the individual components, in this case s, and sy. 
Using the notion that any closed path along the mirror surface must result in a net slope change 


of zero, Freischlad determined the last index of values along the boundaries as [5]: 


N-1 

Sx[N,n] = — y Sx[k,n] (2.4a) 
k=1 
N-1 

sy[m,N] = — )° sy[m,k]. (2.4b) 
k=1 


The two sets of slope data are now NxN. An alternate method is to assume that $[m,n] is 
periodic and that @{1,N +1] = (1, 1] and so on. In this method, the last slopes can be computed, 


also resulting in NxN data. This periodic assumption is included in the definition of the DFT. 


After s, and s, are determined from the sensors, the reconstruction of the b can be completed 
using the FFT2. The first step is to take the FFT2 of the s, and s, data, as 


—1N-1 








ay X silpsale” ) (2.5a) 
p=0 q=0 

N-1N-1 

=y L slp, qle ), (2.5b) 
p=0 q=0 


The data can be filtered in the frequency domain. This filtering is equivalent to solving Equa- 


tions (2.3a) and (2.3b) for @ in the frequency domain to obtain 








A 0, if k,l = 0 
Pk, 1] = jak Bai oes 
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Equation (2.6a) is valid for Hudgin geometry, while Equation (2.6b) is valid for Fried geometry. 


The resulting wavefront phase estimate is available after taking the inverse FFT2: 


1 N-1N-1 


dimnl=— YY Opes m9), (2.7) 


p=0 q=0 


The phase value at the phase points is now known. This information can then be processed by 
the controller to properly deform the mirror surface and improve the resulting image. Since the 
traditional deformable mirror has the actuators co-located with the phase points, the controller 


can easily make the necessary changes that directly improves the imaging results. 


2.2.1 Boundary Conditions 
The combination of the sums in Equations (2.4a) and (2.4b) imply that there are artifacts along 
the boundary of two of the four sides. The sums of noisy data points result in a larger variance 


of the noise on the sum. This is an undesirable consequence of this technique. 


The boundary is where the majority of the error occurs when computing the wavefront. These 
errors can result in large discontinuities which affect the reconstruction algorithm at other points 
across the entire aperture. This can be explained since the impulse response of the algorithm is 
shift variant, although the reconstruction is unbiased [5]. That is to say, the sums in Equations 
(2.4a) and (2.4b) have detrimental effects across the entire aperture in reconstruction; a better 


technique must address this issue. 


2.3 Circular Aperture 

The circular aperture was analyzed by Poyneer [6]. Wavefront reconstruction in the center of 
the aperture is relatively straightforward from the square aperture discussion in Section 2.2. 
However, along the circular boundary, there are large errors in the slopes as calculated by Equa- 
tion (2.3). These occur because the area external to the circular aperture is set to zero in the 
matrix grid of data. The slope data is then calculated incorrectly, and the error propagates into 
the reconstructed wavefront. Several modifications are needed to understand how the slopes are 


affected by the boundary. 


In Figure 2.2, the gradients are shown as scaled vectors sampled at a very low rate. The vectors 
near the boundary contain large errors, as can be seen by the larger magnitude and even opposite 


phase. Additionally, some of the gradients are non-zero outside of the aperture. 


2.3.1 Boundary Conditions 
The definition of the slopes in Equations (2.2) has to be altered to accommodate the circular 


boundary. At each boundary point [m,n], let 


sx|m,n] = ix[m,n| + b,[m,n| (2.8a) 
sy|m,n] = iy[m,n] + by[m,n] (2.8b) 
where the values i,{m,n] and i,[m,n] are internal slopes near the boundary that are measured 
using Equations (2.2); these are the necessary values to compute wave-front phase reconstruc- 


tion correctly. Additionally, the b, [m,n] and b,|m,n] are the boundary slopes that are incorrect; 


removal of these terms is necessary for proper wavefront reconstruction. 


The relationship between the internal and boundary gradients are defined by the matrix 


Mu =c (2.9) 





Figure 2.2: Overhead view of gradients for sample wavefront. 


where M is a matrix determined by the geometry, u is a vector of unknown boundary slopes and 


c is a vector of the measured slopes which slopes cross the boundary. 


The matrix M is developed from the equations that relate the entries of u to c. To solve for u, 
the inverse of M is required. Equation (2.9) has an infinite number of solutions since it does 
not solve for the DC term (colloquially known as the piston). This occurs because, while the 
slopes of the phase are known, the actual phase is not. Since M is singular, M might be better 


conditioned by singular value decomposition while keeping the dominant singular values. 


In Figure 2.3, the unknown boundary slopes are labeled sequentially in u, while the known 
internal slopes are s,{m,n] and sy[m,n]. The u slopes all cross the circular boundary edge and is 
shown as a dotted line. The zero vectors shown indicate that the two connected points are both 


external to the circular aperture. 


By writing the mesh equations for small square edges between adjacent phase points, we obtain 
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Figure 2.3: Circular aperture showing slopes near the boundary. 


a set of linear equations. This method generates the linear combinations necessary to solve for 
the non-noisy data exact solution or the noisy data least-squares solution. These equations di- 
rectly relate the measured slopes to the unknowns and zero external slopes. A starting location 
must be selected along the boundary, and then a path is traced along the boundary until the 
aperture edge is closed. Depending on the geometry, we see that the c entry is a linear combi- 
nation of some s,[m,n| and sy|m,n]. Upon solving for the u solution, we obtain each resulting 
u entry as either a b,{m,n] or by|m,n]. Knowing the b,|m,n] or by[m,n] information, we solve 
Equation (2.8) for the needed i,|m,n] and i,[m,n], which is trivial. The results correct the slope 


data along the boundary, and the data is now used for reconstruction of the phase @. 
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Figure 2.4: SMT wavefront sensor topology for a single segment with modified Fried geometry 
phase point locations. 


2.4 Hexagonal Aperture 

The hexagonal aperture is a modification of a circular aperture. It is more desirable for multiple 
mirror systems than the latter since its segments can easily fit together in a continuous fashion. 
It has been commonly used on many recent space telescopes, including the SMT and the James 
Webb Space Telescope (JWST). 


The linear edges make the mirror segments easier to align with one another, such as by using 
edge-detection sensors. These sensors are able to detect the neighbor segment mirror heights 
and adjust the coarse and fine actuators to place the heights of the mirror surfaces very closely to 
one another at the nanometer scale. This calibration is referred to as phasing the telescope and is 
seldom needed to be repeated for on-orbit operation but is performed often for the ground-based 


telescopes. 


For the SMT, the 61 WES per segment pictured in Figure 2.4 are packed together densely. 
However, their combined edge does not result in an exact hexagonal shape along the outer edge. 
The phase points should be placed along the edges of the individual lenslets. This makes the 
grid uniform in rectangular coordinates. This is a more convenient arrangement than the one 
originally depicted in Figure 2.1. If the phase points were completely inside the lenslet for 
each lenslet, there would be no information about the relationship between the phase points of 
adjacent lenslets. Likewise, if the phase points were all external to the lenslet (as was shown 
in Figure 2.1), there would be an incorrect relationship between phase points and the measured 
slope data due to an overlap between the neighbor lenslet phase points. This arrangement would 


indicate that the phase in one lenslet is influenced by the neighbor, which is not the case. Placing 
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the phase points along the edge in this manner solves these two issues. The geometric impli- 
cations of how these hexagon lenslets focus different amounts of light per axis on the sensor 
detector warrant additional inquiry. The distance between two neighbor phase points across a 
row is different than the distance between two neighbors down a column. This information must 


be included in the phase point calculation, or the visible result will be distorted. 


In the phase point arrangement depicted in Figure 2.4, a single segment of the mirror can be 
stored in a rectangular grid of 10 rows and 18 columns. If all six segments were in a single 


matrix, it can be done in as small of a matrix as 30 rows and 54 columns. 


2.4.1 Boundary Conditions 

The important consequence of the rectangular grid phase point arrangement is that the bound- 
ary correction algorithm for the circular aperture can now be directly applied to the hexagon. 
Without this relationship, the formulas must be modified or the slopes considered to lie along 
the vectors that compose the hexagonal lattice defined by the center points of the 61 WFS. The 
actual formulation of the equations depends on how the aperture boundary crosses the phase 
points. Every possible phase point arrangement should be considered; in a rectangular grid ar- 
rangement of four phase points, there are 16 possible combinations of which phase points are 
inside and which are outside of the aperture boundary; since each of these is handled in exactly 


the same way as they would for the circular aperture, no modification is required. 
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CHAPTER 3: 


Fast Fourier Transform Modeling 





The FFT2 is an implementation of the DFT that is computationally efficient. In a two-dimensional 
grid of @(m,n\, the two-dimensional DFT provides the spatial frequency content of the wave- 
front, showing that the dominant modes of the wavefront are lower frequencies [13]. Therefore, 


in practice only the lower order modes are needed to correct the largest error terms. 


Breaking the wavefront down into spatial modes is well suited for the actuators. On the SMT, 
the actuators are placed into a triangle mesh along the structural supports on the rear surface of 
the mirror. This truss lattice is independent from the hexagon lattice of the WFS. This geometric 
relationship between the phase points and the actuators varies across the entire segment as seen 
in Figures 3.1, C.1 and C.2. There are two methods for dealing with these different lattices. 
The first is a coordinate transformation. The second method is working in the spatial frequency 


domain. 


In the traditional DM controller, the phase points and actuators are co-located. This makes for 
a simple controller because the wavefront phase error is known at the locations of the actuators. 
Thus, a simple actuation change causes a wavefront change, with the goal of being an improve- 
ment in the resulting image. For the SMT geometry, this is not the case. The wavefront phase 


is sampled along the hexagon lattice; the actuators are found along the truss lattice. 
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Figure 3.1: SMT wavefront sensor and actuator topology for a single segment. 
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Figure 3.2: Spring-Mass system model of a deformable mirror. 


3.1 Deformable Mirror Plant Model 


The plant model mathematically represents the physical system, in this case, the DM position. 
To simplify the analysis, the mirror is thought of as a system of masses connected to one another 
with springs. Each mass is connected to four neighbors, as shown in Figure 3.2. In this manner, 
the mirror surface forms a large Cartesian grid of mass points. Each point also has a fifth spring 
connection that is the height displacement out of the plane of the mirror. 


From Figure 3.2, utilizing the concept of Hooke’s law, we can relate the vertical displacement 
y(m,n,t), velocity y(m,n,t) and acceleration j(m,n,t) at location (m,n) and time f to the actu- 
ators u(m,n,t) and the neighbor’s displacement as 


¥(m,n,t) = — ay(m,n,t) — Boy(m,n,t) + You(m,n,t) 
— (0,1(9(m,n,t) —y(m,n—1,1))) 

— (,—1(9(m,n,t) — (m,n + 1,1))) 

— (041,0(9(m,n,t) —y(m—1,n,t))) 

— (Q_1,0(3(m,n,t) —y(m+ 1,n,t))) (3.1) 
— (Bo,1 (v(m, n,t) —y(m,n—1,t))) 

— (Bo,-1(y(m,n,t) — y(m,n + 1,t))) 

— (Bio(y(m,n,t) — y(m—1,n,t))) 

— (B-10(y(m,n,t) —y(m+ 1,n,t))) 


where the & coefficients represent damping and the B coefficients represent the spring effect 
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from the neighbors. The negative signs are to ensure the positive @ and f values result in a 


stable system. 


The transfer function from the input u(m,n,t) to the output y(m,n,t) is defined as 


G(m,n,s) = ae (3.2) 


and the transfer function between the center mass and its neighbors are of the form 


—a(m,n)s — B(m,n) 


se+os+p 2 





H(m,n,s) = 


with m,n = —1,0,1. We also define 6 = %& + 0,1 + Qo,-1 + 1,9 + H1,9 and p = Bo + Bor + 
Bo,-1 + B10 + B_-1,.0. Since o and p are always positive, this system is stable. If @ and B do 
not vary with the index, the result is a symmetric output where the energy dissipates into the 
neighbors evenly. These values affect the pole locations of the denominator, and this determines 
the response time of the system. In this manner, this abstract model can be applied to a DM if 


these values can be determined from either structural analysis or testing. 


Equation (3.1) can be expressed in the time domain as 


1 1 
oy y [A(1,l2,t)y(m — 1 ,n— lh, t —T)] + g(t)u(m,n,t — T) dt (3.4) 


which shows a triple convolution in space and time. In the time-domain analysis, all neighbor 
actuators influence each other. The solution must be obtained by concurrently solving the entire 
set of equations for the entire mirror surface. Equation (3.4) can be written in the Laplace 


domain as 


Y(m,n,s) ey YH (1,,l,8)¥(m—l,n—l,s)| +G(m,n,s)U(m,n,s) (3.5) 
h=-lh=-1 


where G(m,n,s) and H(m,n,s) are given in Equations (3.2) and (3.3). For every s, define the 
two-dimensional NxN FFT2 as 
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—1N-1 





tinge y Y ¥(mn, se “He (km Han) (3.6) 
m=0 n= 
and 
-1N=1 jn 
H(ki, ko, 8) -¥ La H(m,n,s)e 6 (him+kan) (3.7) 
m=0 n= 


with k,,ky =0,...,N—1. Then the discrete spatial convolution in Equation (3.5) becomes 


Y (k1,k2,5) = H(ki,k2,s8)¥ (ki, k2,8) + G(s)U (ki, k2, 8). (3.8) 


This combination of transforms greatly simplifies a three-dimensional convolution into multi- 
plication in the frequency domain. In the frequency domain, each frequency is treated indepen- 
dently from the others; the equations can be solved individually with much less computation 


required. 


3.2 Deformable Mirror Plant Model in State-Space 







Feedforward 
Digits! Filter 
G(k.s) 






Filter 


Feedback 
Digits! Filter 
F(k,s) 


Figure 3.3: Simple block diagram for feed-forward and feedback transfer function. 


The purpose of the mirror controller is to determine the actuator command values at each point 
(m,n) across the entire mirror, u(m,n,t) to remove as much wavefront error as possible. In an 
active-optics system, the actuators might be moved once during a calibration such that the mirror 
form is as accurate as possible, correcting the radius of curvature to match designed radius. In 
an adaptive-optics system, the actuators are used as a function of time to correct for the varying 
wavefront of light. It is possible to attempt to do both of these corrections concurrently, but 


there would be a loss of control authority in that some actuators will have a limited range of 
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motion since the actuators are already using their stroke length for correcting the radius of the 


mirror surface. 


In order to design the controller, first a mathematical model in terms of the state-space equations 
is derived for the simple block diagram of transfer functions shown in Figure 3.3. The general 
relationship between a transfer function and its state-space equations is 


G(s) = =C(s!—A)'B+D (3.9) 





where No(s) is the numerator polynomial and Do(s) is the denominator polynomial of the trans- 
fer function, and A, B, C, D are the matrices that relate the inputs, outputs, and states. In the 


following, the D matrix is removed because it is not used for the model. 


For the feed-forward and feedback blocks, the state-space equations can be derived as 


&, = Az, + dour, (3.10a) 
Vie = COZ, (3.10b) 
5, = AS, + by, (3.10c) 
Wk = CkSj, (3.10d) 
and 
Vk = Ver Wk. (3.11) 


where z, and s, are the internal states of the systems in Figure 3.3. In Equation (3.10), Ao 
relates the neighbor masses to one another, bo relates the actuators to the masses and co relates 
the internal state to the output. Likewise, the coefficients A, b and c, do the same for the 
internal state that varies in spatial frequency. The k subscript indicates that these equations and 
coefficients must be applied to each spatial frequency vector < k,,k2 > individually in actual 
implementation. These equations can be combined to a single state-space model to get 
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x, = Bl (3.12a) 














Sk 
Ao bo 

k= -_ Xp + Uk, 3.12b) 

peg epee | 

Se KS 
=f =s 
(3.12c) 
and 
Yk = leo cs Xj (3.13) 
KS 


=hy 


where f;, g and /,; are the equivalent A, B and C general matrices, respectively. The zero vector, 


Q, indicates its row dimension will match Ag to correctly fill the matrix. 


This state-space model is at the basis of the overall controller we will develop in the next chapter. 
What is significant about this approach is that the models represented by Equations (3.12) and 
(3.13) are all decoupled and can be individually controlled. 
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CHAPTER 4: 


Implementation 





In this chapter, we test the estimation and control techniques developed using simulations in 
MATLAB and Simulink software. Example phase data sets were generated by algorithms de- 
veloped by Melissa Corley [14] while with the Spacecraft Research and Design Center (SRDC) 
at NPS. Her code has excerpts from Chris Wilcox on Spatial Light Modulators (SLM) [15]. 
The simulated phase data size of the aperture was 1024x1024 phase points, much larger than 
existing adaptive-optics systems at the SRDC. This large dimension was chosen to show the 
scalability of the algorithm. 


4.1 Wavefront Reconstruction 
The circular aperture solution forms the basis for solving a number of different geometries. 
Poyneer outlined the general theory of how to accomplish this [6]. In this section, a particular 


algorithm for this technique is developed. 


One of the difficulties of applying a Fourier-based reconstruction algorithm is the fact that 
apertures are generally not square but circular or polygonal, which does not fit with the square 
grid of a FFT2. In order to overcome this difficulty, we need to make some assumptions on the 
data outside the aperture and within the grid of the FFT2. In general, we assume the extra values 
to be zero. Particular care has to be taken at the boundary of the aperture since the approach 


developed by Poyneer leads to algorithms which are fairly sensitive to boundary conditions [6]. 


Apertures with various characteristics of vignetting, such as mechanical obstructions to the op- 
tical pathway, can also be included in the final solution. This can be useful for space-based 
telescope apertures since the secondary mirror usually has mechanical vignetting from its struc- 


tural support that holds the mirror in place. 


To implement this algorithm, the entire grid of data is iterated to identify each group of four 
neighbor phase points. The only groups of consequence are along the boundary. The algorithm 
builds a list of boundary slopes u that needs to be solved. If another new unknown is found 
along the boundary, it is appended sequentially in a list. If it is already in the list, the algorithm 
uses the number previously assigned. For the square, any slopes that are completely internal to 


the aperture need to be used in computing the constant for c. After the entire grid is processed, 
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Figure 4.1: Square of phase points for generating equations. 


all of the equations are identified. The matrix M in Equation (2.9) is quickly produced. All of 
the required information for Equation (2.9) is available to solve the boundary problem. This 


design allows for a variety of apertures to be solved, including the hexagon. 


In Figure 4.1, the black dashed line shows the aperture boundary. The blue phase point dot is 
outside of the aperture, while the black dots are inside. The red dotted line indicates the mesh 
equation direction. The unknowns and internal slopes are identified. In this example, the entry 


in c in Equation (2.9) is equal to sy[j,k + 1] —sy[j+1,k]. 


After correctly setting up Equation (2.9) variables M and c, u is easily solved. The values u 
contains apply to specific s,[j,k] or sy[j,k] that are identified by the unknowns. Each u value is 
iterated and then applied to the slopes in Equation (2.8). The results are matricies of i,| j,k] and 
iy[j,k] which are used in the reconstruction algorithm. This technique is applied to the large 
example data set shown in Figure 1.2. The reconstructed wavefront is shown in Figure 4.2. The 
reconstruction looks comparable to the original phase data, but there are some artifacts in the 


external aperture area. 


Figure 4.2: (a) Original wavefront; (b) Reconstructed wavefront. 
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The difference in the reconstructed wavefront from the original (known) wavefront is computed 
and shown in Figure 4.3. This shows that the larger errors are at the boundary as expected 
or at the external edge of the boundary. This is a consequence of the filtering, discarding the 
imaginary components and round-off errors in floating point operations. Since the algorithm is 
insensitive to the mean value, the comparison is made between the reconstructed data and the 


mean-subtracted original data. 





Figure 4.3: Error in the wavefront reconstruction. 


To better compare the results, the reconstructed data is further processed by zeroing out any 
entries outside of the known aperture. In the actual implementation, these values are not used 
or needed. As can be seen in Figure 4.4, the estimate data closely matches the original data. 
These additional steps are not necessary in actual implementation and are only included here to 


help the reader see that the reconstruction was successful. 


4.2 Fourier Transform Modeling 

The DM is modelled using the discrete spatial Fourier transform (DSFT) to decompose the 
spatial modes of the mirror. Mathematically, this is done using rectangular matrices. Matrix 
sizes are padded with zeros to a length that is a multiple of two. This is done so that the FFT2 
algorithm can execute efficiently. Previous work has been done using 64x64 matrices [6]. In 
this thesis, the matrices used were 256x256 in size to show that the algorithm is scalable to the 


larger apertures that are expected in the future. 


23 


time delay z°. 





Figure 4.4: Cleaned up estimated wavefront reconstruction. 
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Figure 4.5: Plant model for deformable mirror. 


The plant implementation can be seen in Figure 4.5. In order to make the implementation more 
efficient, the analog transfer functions in s are implemented in discrete time in z by converting 
from continuous time to discrete time using a zero-order hold (ZOH) over an assigned sampling 
period 7;. The sampling time can be adjusted to match the physical system response. In order 
to avoid algebraic loops, we implement the transfer functions to be strictly proper by adding a 


1 


The forward gain in Equation (3.2) is implemented as a scalar multiplier of each spatial fre- 


quency. To implement Equation (3.3), two matrices are created, hg, and hg. Each matrix size is 


set to be the same as the size of the aperture and has the form 
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0 — 1 0 0 —Q,-1 
—Q1.0 0 0 0 0 
he(m,n,s) = (4.1) 
0 0 0 0 0 
—A_10 0 0 0 0 


After creating these matrices, their two-dimensional FFTs are computed. The two matrices hg 
and /g are needed for the Laplace s term and the constant term. These matrices are constructed 
with periodicity. The term that relates the mirror mass to its northern neighbor is wrapped 
around to the bottom row of the first column; likewise, the term that relates the mirror mass 
to its western neighbor is wrapped around to the rightmost entry in the first row. The DSFT 
is applied, and the results are identified as matrices Hg and Hg. The three matrices are used 
in conjunction with the digital filters as an element-wise matrix multiplier, which implements 
Equation (3.8) efficiently as compared to Equation (3.4). These matrices are not converted to 


discrete time, as the discrete filters implement the continuous-time equation. 


The purpose of the controller is to meet the design goal of Equation (4.2) given that e(m,n,t) is 
the error of the mirror state minus the sensed wavefront. The ability of the system to keep the 
error below this threshold depends upon how rapidly the mirror can respond to the wavefront’s 


constantly evolving phase progression in time. This requirement is expressed by 


e(m,n,t) = y(m,n,t) — o(m,n,t) (4.2a) 
|e(m,n,t)| < €,Vt >to. (4.2b) 


The last step is to add an integrator state to the final output of the controller. An integrator 
controller has additional robustness compared to classical controllers. The integral controller 
is well suited for tracking of constant reference inputs, which in the case of slowly-varying 


wavefront mode coefficients is an adequate approximation. The resulting state-space equations 


Ux 0 0 Ux 1 = 
= Ee 4.3 
nD allel el <a 
nd 7 
=f, =8 
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and 


ve = [0 hel i] (4.4) 


“sS_ Xk 
= 


For each spatial frequency, the matrix dimensions for f,, @ and hg are 5x5, 5x1 and 1x5, re- 
spectively. However, this fifth order system is not of full rank. An uncontrollable mode exists 
because of a pole/zero cancellation in the transfer function. Referring to Figure 3.3, we obtain 


the transfer function which shows this cancellation: 








Y;(s) a No(s) 1 _ No(s) D(s) a No(s) (4.5) 
Ux(s) ~ Dls) | HG — Ds) DE)—Bils) — DG)— Bus) | 


The pole-zero cancellation is stable and does not affect the stability of the system. The state- 
space equivalent matrices can be computed from the transfer functions which directly define 
the digital filters. The matrices f;, g and h, can all be defined from Equation (4.5) using Equa- 
tion (3.9). Then matrices Fis @ and /y; can be defined from Equations (4.3) and (4.4). At this 
point, the controller is well defined, but an observer needs to be added. The observer estimates 
the mirror height displacement position at each mass position since in the real system the true 
mirror height displacement is not known. The control signal 7 is computed by combining a 
state estimator and state feedback as 


Xp = F Xk + Bie t+ Kelvyn — Xx) (4.6) 


and 


up = —LyX, (4.7) 


with K; such that det(s! — f, + Kyhy) =T(s) and Ly such that det(s/ — f, +2L,) = A(s). Both 
of the polynomials I’(s) and A(s) should be exactly the same for all k in < k,,ko >. In order 
to accomplish this, the values for K;, and Lz must be solved for each k by selecting the pole 
placement locations. These placements can be made to simulate the expected real-world system 
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Figure 4.6: Controller model for deformable mirror. 


with additional information from the SMT. 


All of the system developments up to now lead to the creation of the system shown in Figure 4.6. 
The information needed to produce the combination feed-forward and feedback system that can 
control the DM is highlighted by this model. The observer is made to follow the rule of thumb 
that it converge to the actual state at four to ten times faster than the state changes in the pole 
selection. The number of poles placed in both the controller and observer equals the order 
of the system, which in this case is three. In the figure, the W and V matrices represent the 
polynomial coefficients of the transfer function numerator terms. There are three each, while 


the denominator is a fourth order polynomial that is the same for all k. 


Combining Figures 4.5 and 4.6, we get the resulting complete simulation shown in Figure 4.7. 
The “phi” data input into the system shown in Figure 4.7 is from a phase generation program run 
externally before the simulation [14]. Ideally, in the real system, the wavefront reconstruction 


algorithm discussed in this thesis can be used to generate the estimated phase. 


The simulation is configured with demonstration parameters that exhibit the dynamics expected 
from a DM. The resulting system is stable, and the time-varying wavefront is fed into the system 


controller which provides the actuator displacements to the mirror. 


In Figure 4.8, the first 25 wavefront updates are shown, and the single pixel maximum error is 
plotted. The controller is adjusting the actuators to correct for the error as shown in the plot. The 
wavefront is sampled at a lower rate than the actuators can be moved, so the saw-tooth pattern 
shows that the error is approaching zero. When the new wavefront is provided, the maximum 
pixel error has a discontinuous jump, and the pattern keeps repeating. The percentage is cal- 
culated based on the absolute maximum phase value over the entire aperture. This information 


can be treated as a figure of merit to show that the controller is working across the large set of 
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Figure 4.7: Simulation model for deformable mirror controller. 
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Figure 4.8: Maximum pixel error across the entire aperture for 25 wavefront updates. 


data. 


In Figure 4.9, the pixel error is broken up into the difference between the maximum and the 
minimum. This shows the error envelope. The negative error is much larger initially from the 


large phase changes in the blue regions shown in Figure 1.2. The error on both sides now has 
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Figure 4.9: Maximum error extents on either side of the wavefront across the entire aperture 
for 25 wavefront updates. 


the saw-tooth pattern. At the beginning of the plot, the lower side error makes large jumps from 
larger changes in the example wavefront. Despite the large changes, the system responds by 
a much sharper rate of error decrease. This is an interesting point. The graph results are not 
tracking the same individual pixels throughout the simulation but rather the specific pixels with 


maximum error extents. 


In the spatial frequency domain, the majority of the wavefront energy is concentrated in the 
lower frequencies. This detail can be exploited further to minimize the required computing of 
all spatial frequencies to just the range of interest. If the system can tolerate the higher fre- 
quency content being uncontrolled, the computing requirements will decrease. In this manner, 
extremely large apertures can be processed without higher end computing resources. This can 
also be exploited on a space-based telescope, where the processors lag ground-based perfor- 
mance by nearly a decade. 
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CHAPTER 5: 


Conclusions 





5.1 Summary 

The main focus of this thesis is the application of the SMT to adaptive optics. In particular 
it addressed two important issues: wavefront reconstruction and feedback control. Wavefront 
reconstruction is developed based on the fundamental geometry of the data and the grid points 
where the wavefront phase is measured. Interpolating between these points, we compute the 
estimate of the wavefront. Techniques vary by computational complexity. While ground-based 
telescopes can have large computing power, space-based telescopes usually have other limit- 
ing factors of size, weight, power and electronics reliability to consider. While the industrial 
community is making larger apertures with increased actuator density, the algorithm’s compu- 
tational requirements should not scale faster than the computing capabilities of the platform in 


future systems. 


The algorithms presented in this research are based on the DFT, which is computationally ef- 
ficient and offers results comparable to other techniques such as the Zernike Polynomials. The 


wavefront can be reconstructed for non-noisy data or with a least-squares fit for noisy data. 


Similarly, the control is based on a simple model of the DM as a mesh of masses connected by 
springs. The mirror has actuators along the normal of the rear surface that adjust the heights. 
The model is developed from basic theory and has not yet been adopted yet to the SMT. By 
assuming linearity and stationarity in the time and spatial domains, we can take full advantage 
of a decomposition in terms of spatial frequencies using the FFT2. In this way, the highly 
connected time-space domain model is reduced to a set of independent models in the frequency 


domain which can be easily controlled. 


5.2 Additional Work 


The wavefront reconstruction was accomplished for the single hexagon panel. Additional work 
could apply this to the six panels of the SMT. The interesting contribution of this would be to 
determine if the panels should be solved independently or as a larger more complicated aperture 


to see under which circumstances would the error would improve. 
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In addition, there are a number of tests to run on the actual telescope hardware. The effect of 
mechanical vignetting of the SMT in wavefront reconstruction is not known. The relationship 
of the outer edge of the hexagon aperture to the mirror’s flat edges to the lenslet jagged hexagon 
edges, seen in Figure 3.1, is not known. The lenslets may all be contained inside of the mirror 
edge or extend pass the edge. If there is overlap, it should be studied to determine the effect on 


wavefront reconstruction. 


There has been previous work done by Puschel and Rotteler on the discrete Triangle Trans- 
form (DTT) [16]. This work could be applied to the truss lattice of the actuators. Additionally, 
Vince and Zheng have worked on taking the DFT of a hexagon lattice [17]. This work could be 
applied to the wavefront sensor lattice. Future work could compare which method allows for 
an efficient transformation between lattices and solve the varying relationships for the actuators 
to wavefront sensors in the lattices. There could be interesting consequences of the DTT ap- 
proach given that it is similar to the discrete Cosine Transform (DCT-HI) and is less sensitive to 
boundary conditions. A study of how the data processing requirements vary from this technique 
would be beneficial. 


An interesting consequence of the hexagonal lenslets is that they are not symmetric about the 
two-slope axis. This could have interesting effects on the wavefront reconstruction and sec- 


ondary camera, a topic worth exploring with additional research. 
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APPENDIX A: 
Equation Set for Circular Aperture 





The following equations represent a solution to the boundary conditions represented by Equa- 
tions (2.8) for Figure 2.3: 


ur —u; =0 U3 — Sx[3,2] —u2 =0 

ug + sy[5, 1] — s,[4,2] —u3 =0 us + ug — 5x[5, 2] — sy[5, 1] =0 
5x [6,2] — ug =0 —ug—u7 =0 

ug + ug — $x{7,3] — sy[7,2] =0 —u19 — Ug =0 
U19 — U1 — Sy[8,3] =O Uj, +442 — $x[8,5] — sy[8,4] =0 

5x[8,5] +13 — uy4 — Sy[8,5] = 0 U14 — Uj5 — sy[8,6] = 0 
uy5 —Ujp6 =O 5x[7, 7] +16 — u17 — Sy[7,7] =0 

uj7 —Uujg =O 5x[6, 8] +u1g —uj9 =0 

Sx[5, 8] + uj9 — u29 — Sy[5,8] = 0 5x[4, 8] + sy[5,8] —u21 —u22 =0 
5x[3, 8] + u22 —u23 = 0 ur4 +23 = 0 

5x[2, 7] + sy[3,7] —uo4 — u25 = 0 U6 + U25 = 0 
27 + Sy[2,6] — u26 = O sx[1,5] +.5,[2,5] —u27 — u2g = 0 

u39 + Sy[2,4] — s,[1,5] —u29 = 0 u31 + Sy[2,3] —u39 =0 
—u32 — U3, = 0 uy +$y[3,2] —s,[2,3] —u32 =0 


These are shown for completeness to illustrate the example. These equations are created by 
writing the mesh equation for each square. The path followed along each square is top, right, 
bottom, left. If the slope arrow points in the opposite direction of this path, it is distinguished 


by negative sign. 
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APPENDIX B: 
MATLAB code 





B.1 ComputeSlopes.m 


oe 


This file implements an algorithm for solving the boundary slopes 


ol? 


and then displays the results for a sample set of data. 


oe 


Useful output variables: 


ol? 


unknowns — structure of the unknown slopes. fields are: axis, j, k 


oe 


equations — array of the unknowns for each equation (2 per eqn) 


ol? 


fields are: axis, j, k, sign 


oe 


knowns — array of known slopes needed for each equation (2 per row) 


ae 


This is useful for dynamic, time varying data 


ol? 


Just iterate over this, sum the values per row for new c vector 


ae 


fields are: axis, j, k, sign 


ol? 


M, Minv — matrices that correspond to equations. 


ol? 


code by Travis Axtell 





close all; 











clear all; 





cle; 


% load phi variable 
load data; 


fe) 


% pad zeros on all sides to guarantee an edge 





[sizel, size2] = size(phi); 
phi2 = zeros(sizel+2,size2+2); 
phi2 (2: (sizelt+1l),2:(size2+1)) = phi; 


phi = phi2; 
clear phi2; 


% Find the aperture of data. This code depends on the actual data. 





% For this example data, th dge is entirely surrounded by 0. 
% This should be confirmed before applying new data to algorithm. 


aperture = phi “= 0; 


35 


% Find the slopes from the example phase data. 


% Real world data would just be provided without knowing the phase. 


sx = [phi(:, 2:sizel+2), phi(:, 1)] — phi; 

sy = [phi(2:sizelt+2, :); phi(1l, :)] — phi; 

& Empty variables to begin the code. 

unknowns = []; % list of structures that identify the unknown slopes 
knowns []; % list of structurs that identify the known slopes 

sums = []; % vector of sums, c, to solve for u. 

equations = []; % list of equations using the unknowns 


fe) 


kzero = 


fe) 


for j = 


For 


% Paceholder for zeros as knowns. 


struct. ("“axae", "sero", TU", Jp R*; Ky “Sogn, 1) 7 


% iterate over each square of values from the aperture 


l:sizeltl, 

k = 1l:sizel+l, 

% binary value for the arrangement of the 4 values are as shown: 
1, 2 

4 8 

value = aperture(j,k) + 2x*aperture(j,kt+1) + 


ol? 


ol? 


4xaperture(jt1,k) + 8xaperture(jt1,k+1); 


%& This switch statement is not truly required but helps to organize 
= the code. 

switch value 

case { 0, 15 } 


° 


$ the block is all internal or external and no 





fo) 


% processing is required. 
continue; 
case { 1, 2, 4, 8 } 


fo) 


% single point on aperture. 





% Each if statement below identifies the two unknowns and 
& their associated mesh equation sign. 
sum = 0; 


% kl, k2 are placeholders since they are zeros. 


kl = kzero; 

k2 = kzero; 

if value == 1, 
ul = struct. ("axis*, "x, '7", dy Tht, kK, “ston, 1); 
u2: = struct ("“axas", "yy", "a", dy “et, Kp *sagn'’, —1); 
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elseif value == 2, 








ul = struct ('axis', 

u2 = struct ('axis', 
elseif value == 4, 

ul = struct ('axis', 

uz = struct (“axie", 
elseif value == 8, 

ul = struct ('axis', 

u2 = struct ('axis', 
end 


% checkadd confirms the 


% the unknowns list, or 


is er a "k', Ky "sign', 


47>, jp ‘kt, ke, “sign* 


be air, Jr Paes k, ‘sign’, 


7), SHly RH", ey. *ebgn* 


gt, jt, et) k, Meigat 








14’, jy th", k+l, "sign' 


1); 
, lj); 


—1); 
’ —1); 


’ —1); 
, 1); 


unknown is new before adding it to 


if it is not new, 


redefines the ul 


%& to have the correct unknown number for generating the 


oe 


equations. 


[unknowns, ul] = checkadd(unknowns, ul); 


[unknowns, u2] 


case { 3, 5, 10, 12 } 


% two points on aperture. 


sum = 0; 





checkadd(unknowns, u2); 


% Each if statement below identifies the two unknowns and 





& their associated mesh equation sign. 

if value == 3, 
uh = struct. (“eens Teer, Tat; Ty MMe, Kk, Seagm', 
u2. = struct. ("“axis*, Ty?, '7", gd, ‘ht, kL, “sign' 
sum = sx(j,k); 
kl = struct (*‘axis", "s?, "3", Jy ‘KH, kK, *s2¢en',; 
k2 = kzero; 

elseif value == 5, 
ul. = struct (“exie", "x", "9", dy “Met, Kp Sehan', 
u2 = struct. (taxis; Tt, Va", atl, Th", Kp. ‘sam 
sum = —sy(j,k); % negative for mesh equation 
ki = struct (“eure ,. “ely Tay Jy TH, Ky een", 
k2 = kzero; 

elseif value == 10, 
ul = struct ("axis", "*", "7", Jv Tht, ky, 'Sagn', 
uz = struct.(*axie", x", "7", dtl, "KK", ky “sion 
sum = sy(j,k+1); 
ki = struct (*akia";, "e?, “9% dy, “ER’, kl, “spon 
k2 = kzero; 

elseif value == 12, 


aT 


—1); 
, lj); 


Lye 


1); 


’ -1); 


=1) 7 


116 


117 


118 


119 


120 


121 


122 


123 


124 


125 


126 


127 


128 


129 


130 


131 


132 


133 


134 


135 


136 


137 


138 


139 


140 


141 


142 


143 


144 


145 


146 


147 


148 


149 


150 


151 


152 


153 


154 


155 


156 


ul 
u2 
sum 
k1 
k2 


end 


WP ol? 


oe 


equat 
[unknow 


[unknow 


[o} 


sum = 0; 


° 





& their associated mesh equation sign. 
if value == 7, 
ub. = struct (*exis",; "at, Ta, DHL, Th; Ky “Sign, —1); 
ud = struct.(“axnet,. Nyt, 9", dy “TR, kel, *sugn",. 1) 3 
sum = sx(j,k) — sy(j,k); % signs for mesh equation 
kl = struct (“amas"; "x", "9", dy TE"; Ky *Sagm", 1)¢ 
k2 = struct (*aeag*,. "yt, MA. dy Tet, ky “Segnt) =1); 
elseif value == 11, 
ul -= Struct (“axze", 8, “as gy “ke ">, k, “ehen’,. —1); 
u2. = struct (“axis*, %x?, 7", gl, k"; &k, “szgn*, —1); 
sum = sx(j,k) + sy(j,k+1); % signs for mesh equation 
kil = struct (“axtst, Tat, "7", dy Th'; kp “saign’, 1); 
k? = struct (*exase"™, "yt, “go, ay Me™, Kt, “sign”, 1)¢ 
elseif value == 13, 
ul. = ‘struct: (“axie*, "xs", Ya", dy Te", Ky, “Sngnt, 1) 7 
u2. = struct (texas, “ety. Ty) dy Hh, KEL "Sagm?,;. 1) 4 
sum = —sy(j,k) — sx(jt+1,k); % signs for mesh equation 
ki. =" SExUCE (emia, Tet, Ta ay Tht; Ky “eegm', =1) 7 
k2 = struct (*‘axts", Te?, *7", Jl, "hk", k, *“stgn", —1); 
elseif value == 14, 
ul = struct ("axis*, "*", °3", dy ‘h", ky» “sagn'’;, 1); 
U2. = struct (“axie", “ey, “7, gy “det, k, Seaigiv”, —1); 
sum = —sx(jt+l, k) + sy(j,kt+1); % signs for mesh equation 
kl = struct ("“axis", "x", ' 9’ jel; “e; ky, “sbogn"),. =1)7 


= —sx(jtl, 


= struct ('axis', 


struct (*‘axis", "yy", * 


SEDUCE (tees; THT, 7 


k); %& negat 


ie ane ' 


= kzero; 


the unknowns list, 


j', 
j', 
ive 


7 


i", 


Wy ub Sue 
Vi ret, 


k, 


k+1, 


‘sign’, 


for mesh equation 


j+1, 


or if it is not new, 


WONS*s 
ns, ul] = checkadd(unknowns, ul); 
ns, u2] = checkadd(unknowns, u2); 


T 





case { 7, 11, 13, 14 } 


% three points on aperture. 
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TE, 


k, “sign”, 


—1); 


Meany. 10). 


=1)7 


checkadd confirms the unknown is new before adding it to 


redefines the ul 





to have the correct unknown number for generating the 


% Each if statement below identifies the two unknowns and 


157 


158 


159 


160 


161 


162 


163 


164 


165 


166 


167 


168 


169 


170 


171 


172 


173 


174 


175 


176 


177 


178 


179 


180 


181 


182 


183 


184 


185 


186 


187 


188 


189 


190 


191 


192 


193 


194 


195 


196 


197 


ke? = struct “axas™, "yt, 89%) a> “ie ", KtL, “sign”, 1)¢ 
end 
%& checkadd confirms the unknown is new before adding it to 


% the unknowns list, or if it is not new, redefines the ul 


oe 


to have the correct unknown number for generating the 


% equations. 


[unknowns, ul] checkadd (unknowns, ul); 


[unknowns, u2] checkadd(unknowns, wu2); 


otherwise 


oe 


This code should never get executed since all cases are 





oe 


handled, but is included to verify that the code executes 


oe 


correctly. 


oe 


W: for warning 
fprintf£('W: (4, k) = (%d, %d); value = $%d\n', j, k, value); 


end 


%& The lists are updated to include the new entries 


equations = [equations; ul, u2]; 
sums = [ sums; sum ]; 
knowns = [ knowns; kl, k2 ]; 


end 
end 
% Develop the M matrix from the equations 


egqnien = length (equations) ; 





M = zeros (egnlen,egnlen); 





for m = l:eqnlen, 





% Each row of M contains 1 equation, which includes 2 entries. 


M(m, equations(m,1).number) = equations (m,1).sign; 


M(m, equations (m, 2) .number) equations (m,2).sign; 


end 


& If the Minv is needed to be computed, uncomment this 

SMinv = pinv(M); 

& In this case, the Minv was precomputed. This is how it would be done in 
% actual implementation for faster results. 

load Minv.mat; 


fe) 


% Solving for u is now a simple statement. 
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& For dyanmic data, rebuild sums using knowns before this next line. 


u = Minv*sums; 
% Now have a solution. Just need to apply it to the original data. 
6 sx = ix + bx = ix = sx — bx 





% ix internal slope, sx computed slope, bx boundary slope 


% ix is what is needed, sx is what was measured, bx is the error to remove 


ix = Sx; 
iy = sy; 
unklen = length (unknowns) ; 


for m = l:unklen, 


fo) 


& This if statement determines which set of data to apply the solved 





fo) 


% unknown to. 
if unknowns(m).axis == 'x', 
ix (unknowns (m).j, unknowns(m).k) = 
sx (unknowns (m).j, unknowns(m).k) — u(m); 


elseif unknowns(m).axis == 'y', 


iy (unknowns (m).j, unknowns (m) .k) 
sy (unknowns (m).j, unknowns (m).k) — u(m); 
end 
end 
% The mean is removed from the data before its FFT2 is taken. 
sxfft = fft2(ix—mean(mean(ix))); 


sYfft 


fft2 (iy—mean (mean (iy) )); 


N = sizel+2; 


fo) 


& denom2 is only used to make sure the denom is taking on a good value. 


Sdenom2 = zeros(N,N); 
PHI = zeros(N,N); % preallocated for speed 


for j = 1:N, 
for k = LiN; 





terml = (exp(—1lj * 2 * pi * (k—-1)/N) — 1) * sxXfft(j, k); 
term2 = (exp(—1lj * 2 * pi * (j-1)/N) — 1) * sYfft(j, k); 
denom = 4 « (sin(pi * (j-1)/N)*2 + sin(pi * (k—1)/N)7%2); 


if denom < 0.000000005, 
6 this denominator is getting to be a small value, but we 


% cannot divide by zero. So this checks to prevent that from 


oe 


occurring. Choosing this threshold is based on the N 


ol? 


dimension. For larger N, the threshold must be smaller. 
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239 


240 


241 


242 


243 


244 


245 


246 


247 


248 


249 


250 


251 


252 


253 


254 


255 


256 


257 


258 


259 


260 


261 


262 


263 


264 


265 


266 


267 


268 


269 


fo) 


% This has major changes on the outcome of the 
denom = 0.000000001; 

end 

Sdenom2(j, k) = denom; 


fo) 


%& This PHI equation is from the Poyneer paper. 





reconstruction. 


but 


PHI(j, k) = (terml + term2) / denom; 
end 
end 
% The Poyneer paper specifically states the (1,1) entry needs to be 0, 
& in this filter implementation, the 0 is always there anyways. 
SPHI(1,1)=0; 
% the reconstructed wavefront. Should be completely real but 


*o« some $s 


phinew = 


phimean 





% Error 
figure; 
title ('E 


fo) 


% Origin 


it does have 





mall imaginary components from roundoff errors. 
i1fft2 (PHI); 


= phi-mean (mean (phi) ); 
plot 
imagesc (abs (phimean—phinew) ); 


rror between estimate and actual'); axis off; 


al plot 


figure; imagesc (real (phimean) ); 


title('\ 
% Recons 
figure; 

title('\ 


Phi original data (mean subtracted)'); axis off; 


tructed plot 
imagesc (real (phinew) ) ; 


Phi estimate data'); axis off; 


B.2. CheckAdd.m 


function 
S[list, 
& Select 


% If var 





[ list, entry ] = checkadd( inlist, inentry ) 
entry] = CHECKADD(inlist, inentry) 
ively adds to unknowns list 


iable inentry is contained in the variable inlist, 


& the output is an updated entry with the correct .number 


& and th 





e list is not changed. 
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& If the variable inentry 


fo) 


% new number is assigned. 


len = length(inlist); 
add = 1; 
temp = —1; 


for j = l:len, 
if (inlist(j).axis == 
(inlist(j).j == in 
(inlist(j).k == in 
% this entry alre 
add = 0; 
temp = j; % set t 


end 

end 

if add == 1, 
% add the unknown num 
inentry.number = lent 


fo) 


% set the list to inc 

% the .sign is remove 
list = [inlist, rmfie 

else 

% list is not updated 
list = inlist; 
if temp > -1, 

inentry.number = 

end 

end 

% set entry 

entry = inentry; 


end 


is not in the variable inlist, it is added anda 





The updated list and entry are the output. 


inentry.axis) && 
entry.j) && 
entry.k) 


ady exists, so don't add it. 


he output to the existing entry 


ber 

1; 

lude the new entry 

das it is not important for the list purposes. 


ld(inentry, “sign')]; 


to include new entry 


temp; 
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APPENDIX C: 
Additional Figures 
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Figure C.1: SMT numbered actuator topology for a single segment. 
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